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I.  INTRODUCTION. 


Adjerid  and  Flaherty  [1-3]*  developed  adaptive  finite  element  methods  for  solving 
m-dimensional  vector  systems  of  partial  differential  equations  having  the  form 

d  . 

M(x,/)uf  +  f(x,/,u,Vu)  =  2  [D*(x,/,u)uxJX|,  xefi,  t  >  0,  (la) 

k= 1 

subject  to  the  initial  and  boundary  conditions 

u(x,0)  =  u°(x),  xsQ^dQ,  (lb) 

dm 

either  u{ (x,r )  =  cL(x,t)  or  £  £  u u  (x,r )' vk  =  c, (x,t), 

k=lj=l  ^ 

for  x  €  9Q,  t  >  0,  i  =  1,  2 . m.  (lc) 

They  considered  problems  in  one  (d  =  1)  and  two  (d  =  2)  spatial  dimensions  with 
x  st  [xlf ....  xd]T  denoting  a  position  vector  in  R^,  t  denoting  time,  and  Cl  being  either  a 
segment  of  the  real  line  or  a  rectangle.  The  subscripts  t  and  xk  denote  temporal  and  spa¬ 
tial  partial  derivatives,  respectively,  and  v  =  [V! . vd  ]r  denotes  the  unit  outer  normal 

vector  to  the  boundary  9 Cl  of  Cl.  Problems  were  assumed  to  be  parabolic  and  to  have  an 
isolated  solution;  thus,  M  and  D*,  k  =  1, ...,  d,  are  positive  definite  m  xm  matrices. 

Adjerid  and  Flaherty  discretized  Eq.  (1)  in  space  using  Galerkin’s  method  with  a 
piecewise  linear  polynomial  basis  in  one  dimension  and  piecewise  bilinear  polynomials  in 
two  dimensions.  An  a  posteriori  estimate  of  the  spatial  discretization  error  was  calculated 
using  Galerltin’s  method  with  piecewise  quadratic  functions  in  one  dimension  and  piece- 
wise  cubic  functions  in  two  dimensions.  In  each  case,  a  nodal  superconvergence  property 
of  the  finite  element  method  was  used  to  neglect  errors  at  nodes  and,  thus,  improve  com¬ 
putational  efficiency.  The  error  estimate  was  used  to  control  global  [1]  and  local  [2,  3] 
refinement  procedures  that  added  and/or  deleted  finite  elements  to  the  mesh  in  order  to 

*  References  are  listed  at  the  end  of  this  report. 


satisfy  a  prescribed  global  measure  of  the  spatial  discretization  error.  For  one-dimensional 
problems,  the  error  estimate  was  further  used  to  move  the  finite  element  mesh  so  as  to 
equidistribute  the  global  error  measure.  Ordinary  differential  equations  for  the  finite  ele¬ 
ment  solution,  error  estimate,  and,  in  one  dimension,  mesh  motion  were  integrated  in  time 
using  the  backward  difference  code  DASSL  [18]  for  stiff  differential  and  algebraic  sys¬ 
tems. 

Initially,  a  global  refinement  procedure  was  used  in  combination  with  mesh  motion  to 
satisfy  prescribed  error  tolerances  in  the  Hl  norm  [1].  This  procedure  was  replaced  by  a 
more  efficient  local  mesh  refinement  strategy  and  some  problem  dependent  parameters 
were  removed  from  the  mesh  moving  scheme  [2].  In  particular,  numerical  experiments 
indicated  that  the  performance  of  the  error  estimation  procedure  could  deteriorate  when  the 
system  of  equations  governing  mesh  motion  was  too  stiff.  Adjerid  and  Flaherty  [2] 
remedied  this  defect  by  limiting  the  stiffness  of  the  mesh  moving  equations  and  using 
refinement,  instead  of  mesh  motion,  to  equidistribute  the  error  estimate  in  these  situations. 
They  subsequently  extended  their  finite  element,  error  estimation  and  adaptive  local 
refinement  procedures  to  two-dimensional  parabolic  problems  [3]  and  proved  that  the  error 
estimate  of  References  1  and  2  converged  to  the  true  discretization  error  in  Hl  as  the 
mesh  is  refined  for  linear  one-dimensional  parabolic  systems  [4], 

In  Section  II  of  this  report,  we  review  the  one -dimensional  adaptive  procedures  of 
Adjerid  and  Flaherty  [1,  2],  describe  some  improvements  to  their  mesh  refinement  scheme, 
and  present  some  examples  that  illustrate  the  relationship  and  interaction  between  mesh 
motion  and  refinement  The  essential  details  of  Adjerid  and  Flaherty’s  [3]  two- 
dimensional  procedure  and  the  dynamic  data  structures  used  in  its  implementation  are 
summarized  in  Section  Cl.  The  results  of  a  nonlinear  two-dimensional  example  are  also 
presented  in  Section  ID.  Finally,  in  Section  IV,  we  discuss  our  results  and  suggest  some 
future  directions. 
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II  ONE-DIMENSIONAL  ADAPTIVE  PROCEDURES. 

The  one-dimensional  version  of  Eq.  (1)  consists  of  solving 

MU )uf  +  f(x ,t ,11,11, )  =  [DU ,t ,u)ux ]x ,  x  e  (a, b),  t  >  0,  (2a) 

uU  ,0)  =  u°U),  x  e  [a  M  (2b) 

m 

either  U ,r )  =  c, (t)  or  £  Dij Ujt(x,t)  =  ct (t ), 

M 

forx=a,b,  t  >  0,  i  «  1,  2, ....  m.  (2c) 

The  unit  subscripts  on  x  and  superscripts  on  D  have  been  omitted  for  simplicity. 

The  procedures  for  discretizing  Eq.  (2)  and  estimating  the  spatial  discretization  error 
of  its  solution  are  identical  to  our  earlier  work  [1,  2]  and  are  briefly  summarized  in  Sec¬ 
tion  H.l.  The  essential  details  of  our  current  adaptive  procedure  are  presented  in  Section 

H.2  and  some  examples  illustrating  its  capabilities  and  the  interplay  between  mesh  motion 

and  refinement  are  presented  in  Section  0.3. 

IL1.  Discrete  System.  We  construct  a  weak  form  of  Eq.  (2)  by  assuming  u  6  h£,  select¬ 
ing  a  test  function  v  €  H^,  multiplying  Eq.  (2a)  by  v,  integrating  it  on  a  £x  £  b,  and 
integrating  the  diffusive  term  by  parts  to  obtain 


(v,Mu, )  +  (v,f)  +  A  (v,u)  =  vrDU,;,u)ux  |£,  for  all  v  €  //0l ,  t  >  0, 

where 


(3a) 


b  b 

(v,u)  =  JvU,;)ruU,r)dr,  A  (v,u)  =  JvxrDU,r,u)ux  dx.  (3b,c) 

a  a 

Recall  that  the  Sobolev  space  Hl  consists  of  functions  that  are  square  integrable  and  have 
square  integrable  first  spatial  derivatives.  Functions  belonging  to  Hg  are  further  restricted 
to  satisfy  any  essential  (Dirichlet)  boundary  conditions  in  Eq.  (2c),  while  functions  in  H  0! 
must  satisfy  homogeneous  versions  of  any  essential  boundary  conditions. 
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Initially  u  must  satisfy 

(v,u)  =  (v,u°),  for  all  v  e  H^,  t  =  0,  (3d) 

and  any  natural  (Neumann)  boundary  conditions  in  Eq.  (2c)  should  be  used  to  replace  Dux 
in  the  last  term  of  Eq.  (3a). 

Finite  element  solutions  of  Eq.  (3)  are  constructed  by  selecting  finite  dimensional 
approximations  U  €  Sg  c  Hg  and  Vs  Sq  c  Hq  ofu  and  v,  respectively,  and  finding  U 
such  that 

(V,MUf )  +  (V,f)  +  A  (V,U)  =  \tD(x  ,r,U)Ux  |  *,  for  all  V  e  Sft,  t  >  0, 

(4a) 

(V,U)  =  (V,u°),  for  all  V  6  Sgf,  t  >  0.  (4b) 

Specifically,  we  introduce  a  partition 

7t(r  jV)  :=  [a  =x0(r)  <Xj(r)  <  —  <xN(t)  =  b  }  (5) 

of  [ajb]  into  N  moving  subintervals  (*, _t(r  ),xt-  (r)),  i  =  1,  2 . N,  t  £  0,  and  select  Sjr 

and  S<j  to  consist  of  piecewise  linear  polynomials  with  respect  to  this  partition.  The  sys¬ 
tem  of  ordinary  differential  equations  that  result  from  this  spatial  discretization  can  be 
integrated  in  time  using  one  of  the  many  excellent  software  packages  for  solving  stiff 
differential  systems.  We  found  that  the  backward  difference  code  DASSL  (cf.  Petzold 
[18])  for  differential  and  algebraic  systems  best  fit  our  purposes. 

The  spatial  discretization  error  of  the  finite  element  solution 

e(x,r)  =  u(x,r)  -  U(x,t)  (6) 

satisfies  Eq.  (3)  with  u  replaced  by  U  +  e,  i.e., 

(v,M(U,  +e, ))  +  (v,f(-,r,U+«,U,+ex))  +  A(v.U-t-e)  = 


vrD(;t,/,U+e)(Ux  +  ex)|  *,  for  all  v  €  //<}, 


t  >  0. 


(7a) 


(v,e)  =  (v,u°-U),  for  all  \  e  Hq,  t  =  0.  (7b) 

A  W  All 

We  approximate  e  by  a  function  E  e  S0 ,  where  Sq  is  a  finite  dimensional  subspace  of 
Hq  consisting  of  piecewise  quadratic  functions  that  vanish  on  We  further  approx¬ 

imate  v  by  V  e  Sq  and  determine  E  as  the  solution  of 

(V,M(U,+E,))  +  (V,f(-,r  ,U+E,UX  +EX ))  +  A  (V.U+E)  =  0, 

for  all  V  e  Sq  ,  t  >  0,  (8a) 

(V,E)  =  (V,u°-U),  for  all  V  s  S?,  t  =  0.  (8b) 

In  constructing  the  error  estimate  E(x,r),  we  assumed  the  superconvergence  of  the 
piecewise  linear  finite  element  solution  U(x,r),  i.e.,  we  assumed  that  U(x^)  converges  at  a 
faster  rate  on  7t(r//)  than  elsewhere  on  a  <  x  <  b.  This  superconvergence  property  was 
established  by  Thomee  [19]  and  the  convergence  of  E  to  e  has  been  proven  for  linear 
problems  by  Adjerid  and  Flaherty  [4]. 

The  error  estimate  E(jc  ,/ )  is  used  to  control  the  refinement/coarsening  strategy  and 
the  motion  of  We  determine  mesh  motion  by  solving  the  ordinary  differential  sys¬ 

tem 

i,(0  -  i,-i(D  =  ~UWi  -  W),  i=l,  2 . N,  (9a) 

where  X  is  a  non-negative  parameter,  Wi  is  an  error  indicator  on  the  subinterval  (x,.^,  ), 
and  W  is  the  average  of  Wi%  i  =  1,  2,  ....  N.  We  shall  take  Wi  to  be  the  square  of  the 
local  error  estimate  in  H1,  i.e., 

W,  (r )  =  ||E||,2.,  :=  J  [ErE  +  ExrEx  ]  dx ;  (9b) 

however,  other  local  measures  can  be  used  [2]. 
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When  X  >  0  and  W,-  >  W,  the  right-hand  side  of  Eq.  (9a)  is  negative  and  the  nodes 
X{  and  *,_!  move  closer  to  each  other.  Similarly,  the  nodes  j c,  (r)  and  x^i(:)  move  apart 
when  X  >  0  and  W,-  <  vF.  Coyle  et  al.  [16]  studied  the  stability  of  Eq.  (9a)  with  respect 
to  small  perturbations  from  an  equidistributing  mesh  (i.e.,  one  where  Wi(t)  -  W(t), 
i  =  1,  2, N,  t  £  0)  and  showed  that  such  perturbations  could  only  grow  by  a  bounded 
amount  when  X  >  0,  >  0,  i  =  1,2,  N ,  and  the  velocity  of  the  equidistributing 

mesh  remained  finite  for  t  >  0.  They  further  showed  that  the  mesh  obtained  by  solving 
Eq.  (9a)  stayed  closer  to  the  equidistributing  system,  when  X  was  large.  This,  however, 
introduces  stiffness  into  the  system  which  makes  its  solution  expensive  and,  as  noted, 
causes  seme  difficulties  with  our  error  estimate.  Adjerid  and  Flaherty  [2]  studied  Eq.  (9) 
and  developed  an  adaptive  algorithm  for  selecting  X  as  a  function  of  t  that  balanced 
stiffness  and  equidistribution.  The  procedure  for  selecting  X  will  not  be  discussed  further, 
but  it  has  been  used  in  the  examples  of  Section  H.3. 


In  order  to  maintain  sparsity,  we  eliminate  W  by  combining  Eq.  (9a)  on  two  neigh¬ 
boring  intervals  and  solve  the  scalar  tridiagonal  system 


Xi+i  -  2 xi  +  il+1  =  -\(Wi+i  ~  W'i),  i  =  1,  2 . Af-1,  t  >  0.  (10) 

The  ordinary  differential  equations  resulting  from  Eq.  (8)  and  Eq.  (10)  are  solved 
using  the  same  backward  difference  software  that  is  used  to  integrate  the  finite  element 
system  Eq.  (4). 


IL2.  Adaptive  Algorithms.  In  addition  to  controlling  mesh  motion,  the  error  estimate  E 
described  in  Section  II.  1  is  used  as  an  error  indicator  in  conjunction  with  procedures  that 
locally  refine  or  coarsen  the  mesh.  A  top-level  description  of  an  adaptive  local 
refinement/coarsening  algorithm  is  presented  in  Figure  1  in  a  pseudo-PASCAL  language. 
The  procedure  adapfem  integrates  the  system  Eqs.  (4,  8,  10)  from  time  tinirial  to  [final 
and  attempts  to  keep  the  spatial  error  estimate  ||Ej|i  <  TOL,  where  TOL  is  a  prescribed 
tolerance.  The  time  steps  that  are  selected  by  the  temporal  integration  routine  (e.g.. 


DASSL)  are  denoted  as  Ar[m],  m  =  1,  2, ....  and  the  corresponding  times  are  tout[m], 
m  =  0,  1,  •••,  with  tout  [0]  initially  set  to  tinitial.  The  integration  is  halted  every  nstep 
steps  or  when  tout  [m  ]  =  t final  and  the  arrays  At  and  tout  are  recomputed  with  tout  [0] 
reset  to  the  last  computed  time,  i.e.,  tout  [nstep  ]  or  tfinal . 

The  spatial  error  estimate  ||E||i  is  checked  whenever  the  temporal  integration  is 
halted.  If  HEHj  >  TOL ,  the  last  m  integration  steps  are  rejected  and  the  mesh  is  refined  by 
adding 


N[i]  :=  max  {roundpfljEl^/E]  -  1,  0 }  (11a) 

elements  uniformly  to  (x, •_!,*,),  i  =  1,  2 . N.  Here, 


round  ^(x) 


trunc( x)  +  1,  if  x  -  trunc(x )  £  P 
trunc  (x ),  otherwise , 


where  0  <  P  <  1,  trunc (x )  evaluates  the  integer  part  of  x,  and 


(lib) 


E  :=  Q.9TOL/N.  (11c) 

The  choice  of  P  =  0.2  in  Eq.  (11)  seemed  to  produce  refined  meshes  that  reliably  reduced 
||E||i  to  approximately  TOL  the  next  time  that  the  error  estimate  was  checked.  Further 
justification  for  this  value  of  P  is  given  in  Adjerid  and  Flaherty  [4]. 

The  integration  is  redone  from  rour[0]  to  tout[m  ],  where  m  is  either  nstep  or  such 
that  tout  [m  ]  =  tfinal ,  on  the  refined  mesh  which  has 

Nr=N  +  £N[i]  (12) 

i=i 

elements.  This  process  is  repeated  until  ||E(-,ro«r[m])||i  £  TOL. 


Elements  can  be  deleted  from  a  mesh  whenever  ||E(-,;ou/[m])||i  <  TOL  1 3  or  when¬ 
ever  refinement  was  necessary  to  integrate  from  four[0]  to  tout[m].  The  need  to  refine 


procedure  adapfem  ( tinitial ,  tfinal,  nstep,  TOL ); 

begin 

Calculate  the  initial  conditions  and  an  initial  mesh; 

{  Integrate  the  system  from  tinitial  to  tfinal .  } 

m  :=  0; 

tout  [0]  :=  tinitial ; 
while  tout  [m  ]  <  tfinal  do 
begin 

m  :=  m  +  1; 
redone  :=  false; 

Integrate  Eqs.  (4,  8,  10)  for  one  time  step  At  [m  ]; 
tout [m ]  ;=  tout[m- 1]  +  Ar  [m  ]; 

{  Check  the  error  estimate.  ) 

if  (m  =  nstep)  or  (tout [m ]  =  tfinal)  then 
begin 

Compute  a  new  value  of  X,  if  necessary; 

(  Refine  the  mesh.  } 

while  ||E(  tro«r[m])||1  >  TOL  do 
begin 

Add  elements  to  the  mesh; 

Redo  the  integration  on  the  refined  mesh  from 
t  =  tout  [0]  to  tout[m ]; 
redone  :=  true 
end; 

{  Coarsen  or  regenerate  the  mesh.  } 

if  (||E(-,rour[m])||1  <  TOL  /3)  or  (redone) 
then  Delete  elements  from  the  mesh,  if  possible 
else  Generate  a  new  mesh,  if  necessary; 

tout [0]  :=  tout[m\, 
m  :=  0 

end  {  if  m  =  nstep ...  } 
end  (  while  tout  [m  ]  <  tfinal  } 
end  (  adapfem  }; 


Figure  1.  Top-level  description  of  an  adaptive  finite  element  procedure  with 
mesh  motion  and/or  local  mesh  refinement/coarsening. 


often  indicates  that  the  spatial  error  pattern  has  changed  and  that  fine  grids  may  no  longer 
be  needed  in  some  portions  of  the  domain.  A  mesh  is  coarsened  by  uniting  successive 


pairs  of  elements,  (x,_i,x,)  and  (x,  ,xi+1),  when  ||E(  lrou/[m])||i<y  <  T0L/3N,  j  =  i,  i+1. 
This  union  of  elements  is  only  performed  when  a  significant  percentage  of  elements  may 
be  removed  from  a  mesh.  This  strategy  avoids  the  overhead  associated  with  restarting  the 
temporal  integration  routine. 

If  TOL  /3  £  ||E(-,rouf[m  ])||i  <,TOL,  we  continue  the  temporal  integration  with  the 
existing  mesh  provided  that  its  speed  is  not  too  great  and  it  is  not  close  to  equidistributing 
the  local  error  indicators.  A  mesh  where  the  error  indicators  are  not  equilibrated  indicates 
that  mesh  motion  and/or  refinement  are  being  performed  in  a  suboptimal  manner  and  that 
a  new  mesh  may  be  more  efficient  We  use  the  following  indicator  to  measure  the 
effectiveness  of  a  mesh  n(r  A)  with  respect  to  equidistributing  the  local  error  indicators: 


p(7t(r,N))  =  -j=  £\±Wj-iW\  .  (13) 

WW  istl  jsl 

If  ft  is  a  mesh  that  equidistributes  Wj,  j  -  1,  2, ....  N,  then  Wj  =W,  j  =  1,  2 . N  and 

|i(tc)  =  0.  Larger  values  of  p  indicate  increasing  departures  from  equidistribution.  For 
example,  suppose  all  of  the  error  is  concentrated  in  the  first  element,  i.e.,  W  X=NW  and 
Wj  =  0,  j  =  2,  3, ....  N.  Then  n(?t)  =  N,  which  we  interpret  as  meaning  that  at  least  one 
element  (the  last  one)  will  have  to  cross  N  elements  in  order  to  equidistribute  the  error 
indicators.  If  all  of  the  error  were  concentrated  in  the  N/2  th  element,  then  |i(x)  =  Nil, 


indicating  that  one  element  has  to  cross  N/2  elements  of  n.  Whenever  the  mesh  speed  is 


too  fast  and  p(Jt(rour[/n]//))  >0.1  N,  we  generate  a  new  mesh  that  approximately  equidis¬ 
tributes  the  error  indicators  by  iteratively  removing  elements  with  small  error  indicators 
and  refining  those  having  large  error  indicators. 

Additional  details  of  our  procedures,  such  as  the  generation  of  new  initial  conditions 
whenever  the  number  of  elements  in  a  mesh  changes,  are  as  described  in  Adjerid  and 
Flaherty  [2]. 
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In  Section  H3,  we  present  some  calculations  performed  on  stationary  meshes.  These 
were  done  by  using  a  code  based  on  adapfem  with  the  mesh  moving  parameter  X  =  0. 
Additionally,  we  only  generated  new  meshes  when  ji(7t(tou/[/n]A))  >  0.4 N  in  order  to 
avoid  excessive  restarting  of  the  temporal  integration  routine. 

IL3.  Computational  Examples.  We  conclude  this  section  by  presenting  some  examples 
that  illustrate  our  adaptive  strategies  and  also  attempt  to  appraise  the  relative  advantages  of 
mesh  moving  and  local  refinement  There  are  several  potential  reasons  why  an  adaptive 
procedure  that  combines  mesh  moving  with  refinement  would  be  very  efficient  Mesh 
moving  techniques  are  inexpensive  relative  to  refinement  (cf.  Amey  and  Flaherty  [6])  and 
the  use  of  mesh  motion  should  reduce  the  need  for  refinement  Mesh  motion  can  also 
reduce  the  necessity  of  restarting  the  temporal  integrator,  which  is  an  important  considera¬ 
tion  in  a  method  of  lines  approach  such  as  ours.  Some  refinement  is  essential,  however, 
since  mesh  motion  alone  cannot  generally  satisfy  prescribed  error  tolerances.  Furthermore, 
rapid  mesh  motion,  e.g.,  towards  an  evolving  region  of  high  error,  can  severely  restrict 
time  steps  and  diminish  the  efficiency  of  an  adaptive  procedure  (cf.  Adjerid  and  Flaherty 
[2]).  Finally,  many  more  numerical  techniques  converge  at  higher  rates  on  uniform 
meshes  than  they  do  on  nonuniform  moving  meshes. 

There  is,  thus,  a  need  to  quantify  the  optimal  use  of  mesh  moving  with  local 
refinement;  however,  this  is  a  very  difficult  problem  and  there  have  been  very  few 
attempts  in  this  direction.  Amey  and  Flaherty  [6]  presented  some  computational  results 
comparing  mesh  moving  and  local  refinement  procedures  for  two-dimensional  hyperbolic 
systems.  Bieterman,  Flaherty,  and  Moore  [15]  attempted  to  compare  adaptive  local 
refinement  and  method  of  lines  procedures  for  one-dimensional  parabolic  problems  and 
noted  the  difficulties  in  finding  appropriate  performance  measures.  Herein,  we  apply  a 
code  based  on  our  adaptive  procedures  to  two  computational  examples  and  compare 
results  on  moving  and  stationary  meshes.  We  use  the  total  number  of  space-time  cells  to 
integrate  the  partial  differential  system  from  tinitial  to  t final  as  a  measure  of 


performance.  A  similar  measure  of  computational  complexity  was  used  by  Amey  and 
Flaherty  [6].  It  has  several  apparent  deficiencies,  such  as  not  providing  an  indication  of 
the  effort  devoted  to  the  various  segments  of  the  adaptive  algorithm. 

Example  1.  Consider  the  linear  heat  conduction  problem 


u,  +  u*  +  g(x,t)  =  uxx,  — 1  <  x  <  1,  t>  0.  (14a) 

a(x,0)  =  u°Cx),  -lSxSl,  (14b) 

«(-l,/)  =  C[(/),  u(l,  /)  =  c2(t),  r  2  0.  (14c,d) 

We  select  g,  u°,  clt  Cj,  so  that  the  exact  solution  of  Eq.  (14)  is 

u(x4)  =  1  -  +0-8)1  +  ra/th [20(x  +2r -1.6)]}.  (15) 


Equation  (15)  represents  two  wave  fronts  initially  centered  at  x  =  -0.8  and  x  =  1.6 
and  moving  with  speeds  1  and  -2,  respectively.  The  center  of  the  fastest  front  enters  the 
domain  (-1,1)  at  x  =  1  and  t  =  0.3. 

We  solved  Eq.  (14)  for  0  £  t  £  1.2  using  tolerances  of  2 ~k,  k  =  2,  3,  4,  5,  in  Hl 
with  adaptive  procedures  on  moving  and  stationary  meshes.  The  total  number  of  space- 
time  cells  used  on  0  ^  t  £  1.2,  the  exact  error  \\e  Hi  at  t  -  1.2,  and  the  effectivity  index 

0  -  IIEIMMI1  (16) 

at  t  »  1.2  are  presented  in  Table  1.  The  moving  and  stationary  mesh  trajectories  that  were 
used  to  solve  Eq.  (14)  with  a  tolerance  of  1/8  are  shown  in  Figure  2. 

Solutions  on  moving  meshes  used  less  than  half  of  the  space- time  cells  of  those  on 
stationary  meshes.  A  larger  number  of  cells  are  needed  with  a  stationary  mesh  because 
the  temporal  integration  must  be  restarted  more  often  and  more  time  steps  must  be  redone 
due  to  a  failure  to  satisfy  the  error  tolerance.  In  each  case,  the  actual  error  was  less  than 


the  prescribed  tolerance  and  fine  meshes  were  concentrated  in  high-error  regions.  The 
effectivity  index  is  a  common  method  of  appraising  the  performance  of  an  error  estimation 
technique  (cf.,  e.g.,  Babuska  et  al.  [9]).  Ideally,  0  should  not  deviate  appreciably  from 
unity  and  should  approach  unity  as  N  increases.  The  results  of  Table  1  suggest  that  this 
is  the  case.  The  performance  of  our  error  estimate  seems  to  be  slightly  better  on  a  station¬ 
ary  mesh  than  on  a  uniform  mesh. 


Stationary  Mes 

h 

Movine  Mesh 

Tol. 

No.  Cells 
xKT4 

Ik  111 

e 

No.  Cells 
xlO"4 

Ik  Hi 

6 

1/4 

4.11 

0.1803 

0.983 

1.60 

0.1725 

0.979 

1/8 

8.67 

0.0848 

0.994 

2.79 

0.0903 

0.993 

1/16 

26.94 

0.0413 

0.996 

9.50 

0.0566 

0.988 

1/32 

47.72 

0.0230 

0.998 

20.27 

0.0282 

0.996 

Table  1.  Number  of  space-time  cells  on  0  £  t  £  1.2,  spatial  discretization  error 
at  t  a*  1.2,  and  effectivity  index  at  t  =  1.2  as  functions  of  error  tolerance  using 
stationary  and  moving  mesh  methods  to  solve  Example  1. 


Example  2.  Consider  the  reaction-diffusion  system 

u,  =  Ujgg  -  Du  e-5/r,  LTt  =  +  a Due~*JT , 

0  <  x  <  1,  t  >  0,  (17a,b) 

D  -  Re*laJ5,  (17c) 

u(x,0)  =  r(x,0)  =  1,  0  £  x  £  1,  (17d,e) 

ux(0j)  =  Tx (0,t )  =  0,  ii (1,0  =  7(1,0  =1,  t  >  0.  (17f,g,h,i) 


This  model  was  studied  by  Kapila  [17]  and  used  to  describe  a  single  one-step  reac¬ 
tion  (A  -»  B)  of  a  mixture  in  the  region  0  <  x  <  1.  The  quantity  u  is  the  mass  fraction 
of  the  reactant,  T  is  the  reactant  temperature,  L  is  the  Lewis  number,  a  is  the  heat 
release,  5  is  the  activation  energy,  D  is  the  Damkohler  number,  and  R  >  0.88  is  the  reac- 
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When  L  is  near  unity,  the  temperature  slowly  increases  with  a  "hot  spot"  forming  at 
x  =  0.  At  some  time  /  >  0,  ignition  occurs  and  the  temperature  at  x  =  0  jumps  rapidly 
from  near  unity  to  near  1  +  a.  A  steep  flame  front  then  forms  and  propagates  towards 
x  =  1  with  speed  proportional  to  ga&,2<1+a).  in  practical  problems,  a  is  about  unity  and  8 
is  large;  thus,  the  flame  front  moves  exponentially  fast  after  ignition.  The  solution  tends 
to  a  steady  state  once  the  flame  has  reached  x  -  1. 

We  solved  Eq.  (17)  for  0  £  r  £  0.5  with  a  =  1,  8  =  20,  and  R  =  5  using  tolerances 
of  0.2,  0.1,  and  0.05  on  stationary  and  moving  meshes.  The  number  of  space-time  cells 
needed  to  solve  these  problems  is  presented  in  Table  2  and  the  mesh  trajectories  for  both 
the  stationary  and  moving  mesh  calculations  with  a  tolerance  of  0.1  are  shown  in  Figure  3. 
As  in  Example  1,  the  number  of  stationary  space-time  cells  is  approximately  double  the 
number  of  moving  space-time  cells. 


Table  2.  Number  of  space-time  cells  as  a  function  of  error  tolerance  to  solve 
Example  2  on  0  £  t  £  0.5  using  stationary  and  moving  meshes. 

IH.  TWO-DIMENSIONAL  ADAPTIVE  PROCEDURES. 

Our  finite  element,  error  estimation,  and  local  refinement  procedures  for  two- 
dimensional  partial  differential  systems  closely  parallel  our  one-dimensional  methods  and 
are  briefly  summarized  in  Section  m.l  and  m.2.  The  representation  of  data  and  its 
management  are  much  more  complicated  in  two  dimensions,  and  we  use  a  dynamic  tree- 
data  structure  to  store  and  retrieve  information  about  the  mesh,  solution,  and  error  esti- 


mate.  Similar  structures  have  been  used  by  other  investigators  (cf.,  e.g.,  Babuska  et  al. 
[8-10]  and  Bank  et  al.  [12-14])  to  design  adaptive  procedures  for  elliptic  systems,  and  they 
have  been  shown  to  be  an  effective  means  of  reducing  storage  and  access  overhead.  The 
essential  details  of  our  tree  structure  are  described  in  Section  HI.2  and  a  two-dimensional 
combustion  problem,  similar  to  Example  2,  is  presented  in  Section  m.3. 

HL1.  Discrete  System.  A  weak  form  of  Eq.  (1)  is  constructed  in  the  manner  described  in 
Section  II.  1  for  one-dimensional  problems.  Thus,  we  seek  to  determine  u  e  Hg  such  that 

(v,Uf )  +  (v,f(-,/,u,Vu))  +  A  (v,u)  =  J  [v^D^Vj  +  vrD2ux  vjda, 

da 

for  all  v  e  ,  t  >  0,  (18a) 

(v,u)  =  (v,u°),  for  all  v  e  ,  t  =  0,  (18b) 

where 

(v.u)  =  jf  v(jc  ,y  j  f  u(*  ,y  ,t  )dxldx2,  (18c) 

A  (v,u)  =  J[vJjD1(x,r  ,u)uX|  +  vJjD2(x^,u)uXj]otc  xdx2.  (18d) 

We  have  set  the  mass  matrix  M  in  Eq.  (1)  to  the  identity  matrix  for  simplicity. 

The  functions  u  and  v  are  approximated  by  U  e  Sg  c  Hg  and  V  €  Sq  cHq, 
respectively,  where  Sg  and  Sq  are  spaces  of  bilinear  polynomials  with  respect  to  a  piece- 
wise  rectangular  partition  of  the  rectangular  domain  Q.  The  finite  element  solution  U  is 
obtained  by  solving 

(V.U,)  +  (V,f(-,/,U,VU))  +  A  (V,U)  =  J  [VrD1Uz  V,  + 

da 

for  all  V  €  S  q  ,  t  >0,  (19a) 

(V.U)  =  (V,u°),  for  all  V  €  Sq  ,  t  =  0.  (19b) 


As  in  the  one-dimensional  case,  the  spatial  error  e(x,r)  :=  u(x,r)  -  U(x,r)  is  approxi¬ 
mated  by  E  €  Sf  c///  In  two  dimensions,  we  select  the  finite  dimensional  space  to 
consist  of  piecewise  cubic  functions  with  respect  to  a  piecewise  rectangular  partition  of  Q. 
The  cubic  functions  are  biquadratic  polynomials  that  are  missing  their  quartic  terms  (i.e., 
serendipity  functions  in  the  terminology  of  Zienkiewicz  [20])  and  further  vanish  at  the  ver¬ 
tices  of  each  element.  Thus,  once  again  we  take  advantage  of  nodal  superconvergence  to 
simplify  our  approximation  of  the  discretization  error.  However,  there  is  very  little 
theoretical  justification  of  the  superconvergence  property  for  two-dimensional  problems 
and  we  are  relying  on  computational  evidence  [3]  and  our  one-dimensional  theory  [4]. 

The  approximate  error  E  is  determined  by  replacing  u  and  v  in  Eq.  (18)  by  U  +  E 
and  V €  Sq  c//q  ,  respectively,  where  So  is  composed  of  the  same  cubic  functions  as  Sg, 
and  solving 

(V,  U,+E, )  +  (V ,f(-,r ,U+E,V (U+E)))  +  A  (V,  U+E)  = 

J  [VrD1(U+E)x  vt  +  VrD2(U+E)XlV2]da,  for  all  V  6  Sjf,  t  >  0,  (20a) 

an 

(VJE)  =  (V,u0-U°),  for  all  V  e  sjf,  t  >  0.  (20b) 

The  resulting  ordinary  differential  equations  (19)  and  (20)  for  the  solution  and  error 
estimate  are  integrated  in  time  using  a  code  for  stiff  systems  (e.g.,  DASSL). 

HL2.  Local  Refinement  Algorithms.  A  top-level  description  of  our  two-dimensional 
adaptive  procedure  closely  resembles  the  one-dimensional  algorithm  shown  in  Figure  1, 
except  that  we  have  no  mesh  moving  procedures,  as  yet.  Initially,  the  domain  Q  is  parti¬ 
tioned  into  a  "base"  mesh  of  N  xM  rectangular  elements,  which  is  the  coarsest  mesh  that 
can  be  used  to  solve  the  problem.  Refinement  is  performed  by  bisecting  the  edges  of  a 
coarser  element,  thus,  creating  four  elements  where  there  was  previously  one.  A  base 
mesh  having  four  elements  and  a  refined  mesh  obtained  by  bisecting  one  of  them  is  shown 


in  Figure  4. 


The  refinement  process  may  be  repeated,  i.e.,  elements  may  be  bisected  again  to 
create  four  new  elements.  Additionally,  quartets  of  elements  that  were  created  by 
refinement  may  be  subsequendy  deleted  if  they  are  no  longer  needed  to  maintain  accuracy. 
Bilinear  approximations  in  Sg  and  Sq  and  cubic  approximations  in  Sg  and  Sq  are  con¬ 
strained  to  be  linear  and  quadratic,  respectively,  on  edges  between  elements  of  different 
levels  in  order  to  maintain  continuity  of  U  and  E  on  ft. 

The  mesh  is  organized  as  a  tree  structure  with  the  domain  ft  being  the  root  of  the 
tree  and  the  N  x  M  elements  of  the  base  mesh  being  offsprings  of  the  root  All  nonleaf 
nodes  of  the  tree,  other  than  the  root  node,  have  four  offsprings  which  correspond  to  the 
four  elements  created  by  refining  its  parent  element  The  domain  ft  is  referred  to  as  level 
zero  of  the  tree,  the  elements  of  the  base  mesh  are  level  one,  and  the  levels  increase  as 
elements  are  recursively  refined.  The  tree  structure  for  the  mesh  shown  in  the  lower  por¬ 
tion  of  Figure  4  is  displayed  in  Figure  5. 

Each  node  of  the  tree  contains  the  following  information: 

a.  The  element  number,  say  k ,  of  the  finite  element 

b.  The  level  /  of  the  tree 

c.  Pointers  to  the  four  vertex  nodes  of  element  k 

d.  Pointers  to  the  four  midside  nodes  of  element  k ,  which  are  needed  to  represent  E 

e.  Pointers  to  the  four  elements  neighboring  element  k,  with  a  null  pointer  used  when 
an  edge  of  element  k  is  on  the  boundary 

f.  A  pointer  to  the  parent  of  element  k 

g.  Pointers  to  the  four  sons  of  element  k,  with  null  pointers  used  when  element  k  is  a 


Figure  5.  Tree  representation  of  the  mesh  shown  in  the  lower  portion  of  Figure 


As  in  the  one-dimensional  algorithm  of  Figure  1,  elements  are  added  to  a  mesh  when 
IIEUt  >  TOL  and  deleted  from  a  mesh  when  either  ||E||i  <  TOLt 3  or  when  refinement  was 
necessary  to  integrate  to  the  current  time.  Our  refinement  and  deletion  procedures  impose 
the  following  two  rules,  which  Bank  et  al.  [12-14]  found  to  aid  the  efficiency  and  accu¬ 
racy  of  their  refinement  process  for  elliptic  systems: 

a.  The  1  -irregular  rule,  which  states  that  neighboring  elements  can  differ  by  at  most 
one  level  of  the  tree 

b.  The  3 -neighbor  rule,  which  states  that  any  element  where  the  number  of  edges  con¬ 
taining  elements  at  a  higher  level  of  the  tree  and  the  number  of  boundary  edges  totals 


Refinement  is  performed  by  examining  the  elements  of  a  mesh  by  levels,  proceeding 
from  the  root  to  the  leaf  nodes  of  the  tree.  An  element  k  is  refined  by  dividing  it  into 
four  subelements  whenever  ||E|jijt  >  TOLI^Ne,  where  Ne  is  the  number  of  elements  in 
the  mesh.  Elements  are  deleted  from  meshes,  other  than  the  base  N  xM  mesh,  by  prun¬ 
ing  the  tree.  A  quartet  of  elements  having  the  same  parent  is  deleted  if: 

a.  Every  element  in  the  quartet  has  no  offsprings 

b.  Every  neighbor  of  the  elements  in  the  quartet  are  at  the  same  or  lower  level  of  the 
tree 

c.  The  average  error  estimate  of  the  four  elements  is  less  than  TOL/3yfNj. 

Additional  details  pertaining  to  other  aspects  of  our  adaptive  procedures  are  presented  in 
Adjerid  and  Flaherty  [3]. 

IIL3.  Computational  Example.  A  code  based  on  our  two-dimensional  local  mesh 
refinement  procedure  has  been  written  and  applied  to  several  problems  [3].  Herein,  we 
present  the  results  of  a  two-dimensional  version  of  the  model  combustion  problem  con¬ 
sidered  in  Example  2. 

Example  3.  Consider  the  partial  differential  system  on  the  rectangular  domain 


Q  :=  { I  0  <  xXyx2  <  1  ) 

Tt  =TXiX{  +  TXiX2  +  D(\ +ct-T)e~*n',  (x,y)  e  Q,  t  >  0,  (21a) 

T(x,0)  =  1,  XS  (21b) 

7Xi((U2.0  =  o,  r(U2,r)  =1,  0  £  x2  <,  1,  t  >  0,  (21c,d) 

7,J(x1,0,r)  =  0,  r(x1,l,r)=l,  0  <*!£!,  r  >  0.  (21e,f) 


All  of  the  parameters  are  as  described  in  Example  2.  The  Lewis  number  L  has  been  set 
to  unity  and,  in  this  case,  the  mass  fraction  u  =  1  +  (1  -  T)/ a. 


-22- 


We  solved  Eq.  (21)  with  a  =  1,  8  =  20,  and  R  =  5  using  a  spatial  error  tolerance  of 
0.2.  Mesh  refinement  had  to  be  restricted  to  a  maximum  of  two  levels  because  of  virtual 
memory  restrictions  on  our  computing  system.  The  meshes  that  were  created  at 
t  =  0.2867,  0.2979,  0.3055,  and  1  are  shown  in  Figure  6.  Surface  and  contour  plots  of  the 
calculated  temperatures  at  (  =0.28674  and  0.3115  are  presented  in  Figures  7  and  8, 
respectively. 
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Figure  6.  Meshes  that  were  used  for  Example  3  at  r  =  0.2867  (upper  left), 
0.2979  (upper  right),  0.3055  (lower  left),  and  1  (lower  right). 


The  temperature  slowly  increases  until  ignition  occurs  at  approximately  t  =  0.28. 
The  temperature  at  the  origin  then  jumps  from  near  unity  to  near  two.  A  circularly- shaped 
reaction  front  forms  and  moves  radially  with  a  speed  of  approximately  30  towards  the 
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boundaries.  A  steady  state  is  reached  at  about  r  =  0.32.  Refinement  is  confined  to  the 
vicinity  of  the  reaction  front  The  results  of  Figures  7  and  8  show  some  small  oscillations 
in  the  temperature  ahead  of  the  reaction  front.  At  present,  we  are  unsure  if  these  oscilla¬ 
tions  are  caused  by  interpolation  inaccuracies  in  our  plotting  routines,  inadequate  resolu¬ 
tion  of  the  finite  element  solution  due  to  our  restricting  the  number  of  levels  of  refinement, 
or  an  instability  of  the  reaction  front.  We  plan  to  explore  these  matters  further  using  a 
combination  of  numerical  and  asymptotic  techniques. 

Our  results  on  this  difficult  nonlinear  problem  are  very  encouraging;  however,  we 
anticipate  that  greater  efficiency  could  be  achieved  by  combining  local  mesh  refinement 
with  mesh  moving  as  in  the  one-dimensional  procedures  described  in  Section  II. 

IV.  DISCUSSION  OF  RESULTS  AND  CONCLUSIONS. 

We  have  developed  adaptive  local  mesh  refinement  finite  element  procedures  for 
solving  vector  systems  of  parabolic  partial  differential  equations  in  one  and  two  dimen¬ 
sions.  The  nodal  superconvergence  property  of  the  finite  element  method  on  parabolic 
systems  has  been  used  to  calculate  and  estimate  the  spatial  discretization  error  and  mesh 
motion  has  been  combined  with  local  mesh  refinement  for  one-dimensional  problems. 

Examples  1  and  2  were  designed  to  illustrate  the  performance  of  our  one-dimensional 
procedures  and  to  characterize  the  importance  of  mesh  motion  as  an  adaptive  technique 
relative  to  mesh  refinement.  These  experiments  indicate  that  our  combination  of  mesh 
moving  and  refinement  can  obtain  solutions  with  about  one-half  the  total  number  of 
space-time  cells  of  calculations  performed  using  only  refinement  We  emphasize  the  prel¬ 
iminary  nature  of  these  results.  Many  more  experimental  and  theoretical  investigations 
will  be  necessary  before  firm  conclusions  can  be  reached  regarding  the  optimal  combina¬ 
tion  of  mesh  motion  and  refinement.  Appropriate  performance  measures  and  optimality 
conditions  are  yet  to  be  specified.  There  is  also  a  strong  temptation  to  compare  coded 


implementations  of  procedures  and,  at  this  stage,  we  are  interested  in  more  theoretical 
bounds  on  an  algorithm’s  performance. 

Comparisons  of  the  exact  and  estimated  errors,  presented  in  Example  1  and  in  Refer¬ 
ences  1  through  4,  give  us  some  confidence  in  the  accuracy  of  our  error  estimate.  Addi¬ 
tionally,  the  results  of  Example  3  provide  an  indication  of  the  robustness  of  our  methods. 
This  is  a  difficult  nonlinear  two-dimensional  problem;  yet,  we  solved  it  without  interven¬ 
tion  and  a  priori  knowledge  of  the  solution.  No  special  initial  mesh  was  used,  fine  meshes 
were  automatically  added  to  the  vicinity  of  the  reaction  front,  and  the  fine  meshes  fol¬ 
lowed  the  dynamics  of  the  problem.  Indeed  our  one-  and  two-dimensional  techniques 
seem  to  be  well-suited  for  the  automatic  solution  of  reaction-dL  fusion  systems. 

Despite  our  preliminary  success,  there  is  a  great  deal  more  that  should  be  done  to 
justify  and  improve  the  performance  of  our  procedures.  As  noted  rigorous  analyses  of  the 
convergence  of  our  error  estimate  to  the  true  error  have  only  been  done  for  one¬ 
dimensional  linear  parabolic  problems  on  stationary  meshes  [4].  Dimensional,  nonlinear, 
and  refinement  effects  should  be  included  in  a  complete  analysis.  This  is  a  difficult  task, 
as  very  few  analyses  of  two-dimensional  time  dependent  problems  with  refinement  have 
appeared  in  the  literature. 

Several  computational  procedures  in  our  approach  might  also  be  improved.  For 
example,  a  sparse  Gaussian  elimination  procedure  was  used  to  solve  the  linear  algebraic 
systems  associated  with  the  temporal  integration  of  Eqs.  (4),  (8),  and  (10)  in  one  dimen¬ 
sion  and  Eqs.  (19)  and  (20)  in  two  dimensions.  The  solution  of  linear  systems  is  a 
significant  pan  of  the  total  computational  effort,  and  it  is  possible  that  iterative  schemes, 
such  as  multigrid  methods,  could  substantially  improve  performance  and  reduce  storage. 
Multigrid  iteration  was  used  successfully  in  the  adaptive  PLTMG  package  for  elliptic  sys¬ 
tems  by  Bank  et  al.  [13]. 


We  are  also  studying  the  addition  of  mesh  moving  capabilities  to  our  two- 
dimensional  algorithm,  the  use  of  higher-order  finite  element  approximations,  and  imple¬ 
mentations  of  our  procedures  on  vector  and  parallel  computers.  A  simple,  stable  and 
explicit  mesh  moving  technique,  that  may  be  useful  for  our  purposes,  was  developed  by 
Amey  and  Flaherty  [5]  for  two-dimensional  hyperbolic  systems.  This  procedure  dramati¬ 
cally  reduced  errors  and  enhanced  the  resolution  of  their  solutions  (cf.  Amey  and  Flaherty 
[6]).  We  are  developing  procedures  for  two-dimensional  parabolic  problems  that  use 
piecewise  biquadratic  finite  element  approximations  as  solution  spaces  and  piecewise  cubic 
approximations  as  error  estimates.  Babuska  [7]  has  shown  that  the  error  associated  with 
even-degree  polynomial  finite  element  approximations  for  elliptic  problems  is  principally 
due  to  the  error  in  the  interior  of  the  element  Thus,  the  error  on  element  boundaries  may 
be  neglected.  Babuska  and  Yu  [11]  have  implemented  procedures  for  elliptic  systems 
based  on  this  theory  and  we  are  studying  their  utility  for  parabolic  problems.  Finally,  our 
tree  structure  is  well-suited  for  parallel  computation  and  we  are  exploring  its  use  on  a 
variety  of  parallel  computing  systems. 


-28- 


REFERENCES 

1.  S.  Adjerid  and  J.E.  Flaherty,  “A  Moving  Finite  Element  Method  with  Error  Estima¬ 
tion  and  Refinement  for  One-Dimensional  Time  Dependent  Partial  Differential  Equa¬ 
tions,”  SIAM  J.  Numer.  Anal.,  23  (1986),  pp.  778-796. 

2.  S.  Adjerid  and  J.E.  Flaherty,  “A  Moving  Mesh  Finite  Element  Method  with  Local 
Refinement  for  Parabolic  Partial  Differential  Equations,”  Comp.  Meths.  Appl.  Mech. 
Engr.,  56  (1986),  pp.  3-26. 

3.  S.  Adjerid  and  J.E.  Flaherty,  “A  Local  Refinement  Finite  Element  Method  for 
Two-Dimensional  Parabolic  Systems,”  Tech.  Rep.  No.  86-7,  Department  of  Com¬ 
puter  Science,  Rensselaer  Polytechnic  Institute,  Troy,  1986. 

4.  S.  Adjerid  and  J.E.  Flaherty,  “Local  Refinement  Finite  Element  Methods  on  Station¬ 
ary  and  Moving  Meshes  for  One-Dimensional  Parabolic  Systems,”  in  preparation. 

5.  D.C.  Amey  and  J.E.  Flaherty,  “A  Two-Dimensional  Mesh  Moving  Technique  for 
Time  Dependent  Partial  Differential  Equations,”  Tech.  Rep.  No.  85-9,  Department  of 
Computer  Science,  Rensselaer  Polytechnic  Institute,  Troy,  1985.  Also  J.  Compux. 
Phys.,  67  (1986),  pp.  124-144. 

6.  D.C.  Amey  and  J.E.  Flaherty,  “An  Adaptive  Method  with  Mesh  Moving  and 
Refinement  for  Time-Dependent  Partial  Differential  Equations,”  Trans.  Fourth  Army 
Conf.  Appl.  Maths,  and  Comput.,  ARO  Report  87-1,  U.  S.  Army  Research  Office, 
Research  Triangle  Park,  NC,  (1987),  pp.  1115-1142. 

7.  I.  Babuska,  personal  communication,  1986. 

8.  I.  Babuska  and  A.  Miller,  “A  Posteriori  Error  Estimates  and  Adaptive  Techniques  for 
the  Finite  Element  Method,”  Tech.  Note  BN-968,  Institute  for  Physical  Science  and 
Technology,  University  of  Maryland,  College  Park,  MD,  1981. 


■*!  ■  A  »'t 


WWU WW1I 


Mmuwii»ujnwu*i  *\iwMwvw\wwwvwwv  wmw^wwi. 


-  29- 

9.  I.  Babuska,  A.  Miller,  and  M.  Vogelius,  “Adaptive  Methods  and  Error  Estimation  for 
Elliptic  Problems  of  Structural  Mechanics,”  in  Adaptive  Computational  Methods  for 
Partial  Differential  Equations ,  I.  Babuska,  J.  Chandra,  J.E.  Flaherty,  eds.,  SIAM,  Phi¬ 
ladelphia,  1983,  pp.  57-73. 

10.  I.  Babuska  and  W.  Rheinboldt,  “Error  Estimates  for  Adaptive  Finite  Element  Com¬ 
putations,”  SIAM  J.  Numer.  Anal.,  15  (1978),  pp.  736-734. 

11.  I.  Babuska  and  D.  Yu,  “Asymptotically  Exact  A-Posteriori  Error  Estimator  for 
Biquadratic  Elements,”  Tech.  Note  BN-1050,  Institute  for  Physical  Science  and 
Technology,  University  of  Maryland,  College  Park,  MD,  1986. 

12.  R.E.  Bank,  “PLTMG  Users’  Guide,”  June  1981  version.  Technical  Report,  Depart¬ 
ment  of  Mathematics,  University  of  California  at  San  Diego,  La  Jolla,  CA,  1982. 

13.  R.E.  Bank  and  A.  Sherman,  “An  Adaptive  Multi-Level  Method  for  Elliptic  Boundary 
Value  Problems,”  Computing,  26  (1981),  pp.  91-105. 

14.  R.E.  Bank,  A.H.  Sherman,  and  A.  Weiser,  “Refinement  Algorithms  and  Data  Struc¬ 
tures  for  Regular  Local  Mesh  Refinement,”  Mathematics  and  Computers  in  Simula¬ 
tion,  to  appear. 

15.  M.  Bieterman,  J.E.  Flaherty,  and  P.K.  Moore,  “Adaptive  Refinement  Methods  for 
Non-Linear  Parabolic  Partial  Differential  Equations,”  Chap.  19  in  Accuracy  Esti¬ 
mates  and  Adaptive  Refinements  in  Finite  Element  Computations,  I.  Babuska,  O.C. 
Zienkiewicz,  J.R.  Gago,  and  E.R.  de  A.  Olivera,  Eds.,  John  Wiley  and  Sons,  Chiches¬ 
ter,  1986. 

16.  J.M.  Coyle,  J.E.  Flaherty,  and  R.  Ludwig,  “On  the  Stability  of  Mesh  Equidistribution 
Strategies  for  Time-Dependent  Partial  Differential  Equations,”  J.  Comput.  Phys.,  62 
(1986),  pp.  26-39. 


,1 


-  30- 


17.  A.K.  Kapila,  Asymptotic  Treatment  of  Chemically  Reacting  Systems ,  Pitman 
Advanced  Publishing  Program,  Boston,  1983. 

18.  L.R.  Petzold,  “A  Description  of  DASSL:  A  Differential/ Algebraic  System  Solver,” 
Sandia  Report  No.  Sand.  82-8637,  Sandia  National  Laboratory,  Livermore,  CA,  1982. 

19.  V.  Thomee  ,  ‘‘Negative  Norm  Estimates  and  Superconvergence  in  Galerkin  Methods 
for  Parabolic  Problems,”  Math.  Comp.,  34  (1980),  pp.  93-113. 

20.  O.C.  Zienkiewicz,  The  Finite  Element  Method:  Third  Edition,  McGraw  Hill,  London, 


TO 


lUfj 


ww 


WWWffTOWWI WWW?  fWWWW  W  M  W!  WTWWTWW  yrgf  wn. 


TMT-HTT 


'"cCHNICAL  REPORT  INTERNAL  DISTRIBUTION  LIST 


CHIEF,  DEVELOPMENT  ENGINEERING  BRANCH 
ATTN:  SMCAR-CC3-D 
-DA 
-DC 
-DM 
-DP 
-DR 

-DS  (SYSTEMS) 


CHIEF,  ENGINEERING  SUP»OR"  3PANCH 
ATTN:  SMCAR-CC8-S 

-SE 

CHIEF,  RESEARCH  BRANCH 
ATTN :  SMCAR-CCB-R 

-R  (ELLEN  FOGARTY) 

-RA 

-RM 

-RP 

-RT 


i 

r- 

r- 

r- 

t- 

! 


TECHNICAL  LIBRARY 

ATTN:  SMCAR-CCB-TL 

TECHNICAL  PUBLICATIONS  &  EDITING  UNIT 
ATTN:  .  SMCAR-CCB-TL 

DIRECTOR,  OPERATIONS  DIRECTORATE 
ATTN :  SMCWV-00 

DIRECTOR,  PROCUREMENT  DIRECTORATE 
ATTN:  SMCWV-PP 

DIRECTOR,  PRODUCT  ASSURANCE  DIRECTORATE 
ATTN:  SMCWV-QA 


NOTE: 


NO.  CF 
COPIES 


1 


1 

1 

1 


1 


2 

1 

1 

1 

1 

5 


2 


1 


PLEASE  NOTIFY  DIRECTOR,  BENET  WEAPONS  LABORATORY,  ATTN:  SMCAR-CC3-TL , 
OF  ANY  ADDRESS  CHANGES. 


TOBEnropnBTipnpnen^ticTwnongnPTWTignwnignpncng^^ 


TECHNICAL  REPORT  EXTERNAL  DISTRIBUTION  LIST 


NO.  OF 
COPIES 


NO.  OF 
COPIES 


ASST  SEC  OF  THE  ARMY 
RESEARCH  AND  DEVELOPMENT 
ATTN:  DEPT  FOR  SCI  AND  TECH 
THE  PENTAGON 

WASHINGTON,  D.C.  20310-0103 
ADMINISTRATOR 

DEFENSE  TECHNICAL  INFO  CENTER 
ATTN:  OTIC-FDAC 
CAMERON  STATION 
ALEXANDRIA,  VA  22304-6145 

COMMANDER 
US  ARMY  ARDEC 
ATTN:  SMCAR-AEE 

SMCAR-AES,  BLDG.  321 
SMCAR-AET-O,  BLDG.  351N 
SMCAR-CC 
SMCAR-CCP-A 
SMCAR-FSA 
SMCAR-FSM-E 
SMCAR-FSS-D,  BLDG.  94 
SMCAR-MSI  (STINFO) 
P1CATINNY  ARSENAL.  NJ  07806-5000 


12 


COMMANDER 

ROCK  ISLAND  ARSENAL 

ATTN:  SMCRI-ENM 

ROCK  ISLAND,  IL  61299-5000 

DIRECTOR 

US  ARMY  INDUSTRIAL  BASE  ENGR  ACTV 
ATTN:  AMXIB-P 

ROCK  ISLAND,  IL  61299-7260 
COMMANDER 

US  ARMY  TANK-AUTMV  R&D  COMMAND 
ATTN:  AMSTA-DDL  (TECH  LIB) 

WARREN,  MI  48397-5000 

COMMANDER 

US  MILITARY  ACADEMY 

ATTN:  DEPARTMENT  OF  MECHANICS 

WEST  POINT,  NY  10996-1792 

US  ARMY  MISSILE  COMMAND 
REDSTONE  SCIENTIFIC  INFO  CTR 
ATTN:  DOCUMENTS  SECT,  BLDG.  4484 
REDSTONE  ARSENAL,  AL  35898-5241 


ctr 


DIRECTOR 

US  ARMY  BALLISTIC  RESEARCH  LABORATORY 
ATTN:  SLCBR-DD-T,  BLDG.  305  1 

ABERDEEN  PROVING  GROUND,  MD  21005-5066 

DIRECTOR 

US  ARMY  MATERIEL  SYSTEMS  ANALYSIS  ACTV 
ATTN:  AMXSY-MP  1 

ABERDEEN  PROVING  GROUND,  MD  21005-5071 

COMMANDER 
HQ,  AMCCOM 

ATTN:  AMSMC-IMP-L  1 

ROCK  ISLAND,  IL  61299-6000 


NOTE:  PLEASE  NOTIFY  COMMANDER,  ARMAMENT  RESEARCH,  DEVELOPMENT,  AND  ENGINEERING 
CENTER,  US  ARMY  AMCCOM,  ATTN:  BENET  WEAPONS  LABORATORY,  SMCAR-CCB-TL , 
WATERVLIET,  NY  12189-4050,  OF  ANY  ADDRESS  CHANGES. 


COMMANDER 

US  ARMY  FGN  SCIENCE  AND  TECH 
ATTN:  DRXST-SD 
220  7TH  STREET,  N.E. 
CHARLOTTESVILLE,  VA  22901 

COMMANDER 
US  ARMY  LABCOM 
MATERIALS  TECHNOLOGY  LAB 
ATTN:  SLCMT-IML  (TECH  LIB) 
WATERTOWN,  MA  02172-0001 


TECHNICAL  REPORT  EXTERNAL  DISTRIBUTION  LIST  (CONT'D) 


NO.  OF 
COPIES 

COMMANDER 

US  ARMY  LABCOM,  ISA 

ATTN:  SLCIS-IM-TL  1 

2800  POWDER  MILL  ROAD 
AOELPHI,  MD  20783-1145 

COMMANOER 

US  ARMY  RESEARCH  OFFICE 

ATTN:  CHIEF,  IPO  1 

P.O.  BOX  12211 

RESEARCH  TRIANGLE  PARK,  NC  27709-2211 
DIRECTOR 

US  NAVAL  RESEARCH  LA8 

ATTN:  MATERIALS  SCI  &  TECH  OIVISICN  1 

CODE  26-27  (DOC  LIB)  1 

WASHINGTON,  D.C.  20375 


NO.  OF 
COPIES 

COMMANDER 

AIR  FORCE  ARMAMENT  LABORATORY 

ATTN :  AFATL/MN  1 

EGLIN  AFB ,  FL  32543-5434 

COMMANDER 

AIR  FORCE  ARMAMENT  LABORATORY 

ATTN:  AFATL/MNG 

EGLIN  AFB,  FL  32542-5000 

METALS  AND  CERAMICS  INFO  CTR 
BATTELLE  COLUMBUS  DIVISION 
505  KING  AVENUE 
COLUMBUS,  OH  43201-2593 


NOTE PLEASE  NOTIFY  COMMANDER,  ARMAMENT  RESEARCH,  DEVELOPMENT,  AND  ENGINEERING 
CENTER,  US  ARMY  AMCCOM,  ATTN:  BENET  WEAPONS  LABORATORY,  SMCAR-CCS-TL , 
WATERVLIET,  NY  12189-4050,  OF  ANY  ADDRESS  CHANGES. 


